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Bulk fermionic matter, as it can be notably found in supernova matter and neutrons stars, is 
subject to correlations of infinite range due to the antisymmetrisation of the N-body wave function, 
which cannot be explicitly accounted for in a practical simulation. This problem is usually addressed 
in condensed matter physics by means of the so-called Twist Averaged Boundary Condition method. 
A different ansatz based on the localized Wannier representation has been proposed in the context 
of antisymmetrized molecular dynamics. In this paper we work out the formal relation between the 
two approaches. We show that, while the two coincide when working with exact eigenstates of the N- 
body Hamiltonian, differences appear in the case of variational approaches, which are currently used 
for the description of stellar matter. Some model applications with Fermionic Molecular Dynamics 
are shown. 

PACS numbers: 21.60.-n 26.60.-c 71.10.Ca 



I. INTRODUCTION 



Interacting fermionic systems in the bulk limit are a standard object of theoretical study in condensed matter 
physics. Electrons are subject to an external periodic potential in the presence of a cristalline ionic structure and 
the bulk limit can be seen as an infinite number of spatial replicas of a finite system within a specific geometry [l|. 
In that case, the observables of the bulk system can be obtained from the modelization of one single elementary cell, 
provided adequate boundary conditions are applied to the many-body wave function. This amounts to introduce a 
Bloch phase or twist to each wave function in the single-particle basis, and average the twisted observables over the 
different phases within the first Brillouin zone0, In practical applications, this technique has been applied to 
Quantum Monte-Carlo simulations of electron systems also in the absence of any external periodic potential[3-13l- 
this case, tlie introduction of Bloch phases has to be understood as a technique to accelerate the convergence towards 
the thermodynamic limit of these very expensive numerical calculations, which would otherwise become prohibitive 
in computation time. In the absence of an external potential, the periodic cell is just a computation cell with no 
physical meaning, and independence of the results respect to its size has to be checked. 

Coming to the strongly interacting fermionic systems studied in nuclear physics, the bulk limit has not attracted 
much interest in the community since nuclei are finite. However there are physical situations where the propertis 
of fermionic matter composed of protons and neutrons in the thermodynamic limit are essential such as for core- 
collapsing supernova, and for the crust of the neutron stars which are left over by the explosion. This nucleonic stellar 
matter covers a very wide domain of densities ranging from p k, 10* g • cm~^ to a few times the normal saturation 
nuclear density p ~ 10^^ g • cm^'^, temperatures between less than 1 and more than 20 MeV, and proton fractions 
varying between 0. and 0.5 In the sub-saturation density regime, it is well established that matter is charge neutral 
and mainly composed of neutrons, protons, electrons, positrons and photons in thermal and typically also chemical 
equilibrium Q . Depending on the thermodynamic condition, neutrinos and anti- neutrinos can also participate to 
the equilibrium. 

A very large amount of literature exists on the microscopic modelizations of the neutron star outer and inner 
crust |7H12|. In this regime at zero temperature and relatively low density, matter consists of a lattice of Wigner-Seitz 
cells, each cell containing a spherical neutron-rich nucleus immersed in a sea of dilute gas of neutrons and relativistic 
electrons uniformly distributed inside the cell[l^ [l^. The linear size of the cell is of several hundreds of fermis in 
the outer crust. The size decreases going towards the center of the neutron star and the central nucleus becomes 
heavier and increasingly neutron rich. Typical values[13] in the inner crust range from about 200 particles in a cell 
of linear size of about 50 fm for p k, 10^^/m^'^ to about 1500 particles in a cell of linear size of about 15 fm for 
p « 8 • lO^^/m""^. Most of the existing calculations are based on variational approaches (Hartree-Fock or Hartree- 
Fock-Bogoliubov) and employ mixed Dirichlet-Neumann boundary conditions at the edge of the cellfl^, chosen to 
produce a flat density close to the cell border and thus to simulate a uniform neutron gas. The specific way of fixing 
these mixed boundary conditions is not completely clear, and discrepancies due to the choice of fixing these conditions 
increase with density [15.]. A more conceptual problem with Dirichlet or Neumann boundary conditions, is that they 
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neglect antisymmetrization correlations which extend beyond the cell size. This has been pointed out in refs.[16|-[l9j, 
where Bloch boundary conditions, similar to the ones used for the QMC modelizations of bulk electron systems, have 
been employed. In these works it was shown that both the neutron specific heat and the motion of unbound neutrons 
are affected by these effects. 

In the intermediate region between crust and core, complex phases are predicted that can break the spherical 
symmetry of the Wigner-Seitz cell, and can violate the associated translational invariance. This is even more true at 
finite temperature, where the periodicity of the Wigner-Seitz cell is broken by thermal agitation. These conditions are 
met in stellar matter in the pre- and post-bounce supernova dynamics, as well as in the cooling process of proto-neutron 
stars. Even when periodicity cannot rigourously be assumed in these thermodynamic conditions, it can be kept as a 
practical working hypothesis allowing to address the thermodynamic limit. As already mentioned, the drawback of 
that is that convergence with respect to the cell size has to be systematically checked. Time-dependent variational 
microscopic calculations have bee n prop osed to address this region, where statistical averages are calculated from 
time averages assuming ergodicitv |20l - [23| . These works have shown that structures, though they can be degenerate 
in energy [23j, are nonetheless approximately periodic in space even at finite temperature. 

In these calculations simple periodic boundary conditions are employed, completely neglecting the antisymmetriza- 
tion correlations beyond the calculation grid. The importance of properly accounting for these correlations was 
recently stressed in ref. 24] . In this work, a specific ansatz for the boundary conditions based on the localized Wannier 
representation has been proposed in the context of antisymmetrized molecular dynamics. It was shown that such 
boundary conditions allow obtaining a distribution similar to the one of a free Fermi gas, if gaussian wave packets of 
fixed width are periodically disposed on a two-dimensional grid, while an artificial Pauli potential is needed in order 
to obtain the same result with classical molecular dynamics. This result implies that, in the inner crust region where 
stellar matter contains an important component of quasi-free neutrons, properly accounting for the periodic character 
of the system may be of importance. 

To conclude these introductory remarks, it appears that simple Dirichlet, Neumann, or periodic boundary conditions 
are not adapted to the description of bulk fermionic matter. Different solutions are proposed in the literature for 
different applications, naniely the Twist Averaged Boundary Conditions (TABC) in QMC0,|1], the Bloch method 
in mean- field calculations 1 1 6l - [l9| . and the Wannier replica method for molecular dynamics approaches (24j. but the 
equivalence of the different techniques and/or their domain of validity is not completely clear. 

In this paper, we will formally develop the link between the different methods and show some model applications 
for simple non-interacting nuclear systems described through the variational Fermionic molecular Dynamic (FMD) 
method. We will show that Bloch (or TABC) and the replica method are equivalent when they are applied to the exact 
eigenstates of the many-body Hamiltonian. When this is not the case, as for variational mean-field theories applied 
to interacting systems, the equivalence is broken. In this case, the Bloch or TABC technique can produce solutions 
which differ from the exact results more than if simple periodic boundary conditions are applied. Conversely, the 
replica method appears more powerful and can easily be applied to any variational based mean-field treatment with 
an affordable extra computational cost. Our numerical applications will concern one-dimensional systems, for which 
it has been recently argued[25'] that additional drawbacks appear in the TABC method. However it is important to 
remark that one-dimensional modelizations can often be used in the stellar matter case where (except at very high 
density close to saturation) spherical symmetry is often a good approximation. 

II. THE BLOCH THEOREM AND TWIST AVERAGED BOUNDARY CONDITIONS 

We want to address the physical problem of an infinite system (specifically; the neutron star crust) constituted of 
an infinite number of spatial replicas of finite systems of linear size L (the Wigner-Seitz cells). In the absence of any 
periodic external potential, for the replicated system to be equivalent to the infinite system, each Wigner-Seitz cell 
is supposed to contain an integer number of structures (spherical nuclei or more exotic 'pasta' structures [20l - [23| ) of 
size I J Li — Tiili^ i ^ x,y, z. The minimal choice is to take a box containing exactly one structure, and in this case 
the symmetry of the box will follow the symmety of the physical system (a cylindrical box for cylindrical structures, 
etc). To simplify the notations we will work in one dimension only and write L = nl. The extension to three dimensions 
is straightforward. 

A. The Bloch theorem at the N-body level 

We are interested in the translational invariance properties of the Hamiltonian induced by the imposed periodicity. 
The Hamiltonian of the global system is obviously invariant respect to a simultaneous translation of all particle 
coordinates of the same arbitrary length r. This invariance physically corresponds to the general statement that the 
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center of mass momentum is a good quantum number. This symmetry has no influence in an infinite system, and we 
will not consider it further. The fact of working in a finite box of linear size L which is replicated induces additionally 
an extra non trivial invariance, which we now discuss. The total Hamiltonian of the infinitely replicated system reads 



H=Y.Hl (^(™)) (1) 



where x^™^ = (x™, . . . , x^) denotes coordinates of particles belonging to the m-th replica and the cell Hamiltonian is 

N oo N 

i—1 m— — oo i>j 

Note that the cell Hamiltonian depends on a finite number N of particles even if these particles are not necessarily 
confined into the specific volume of the cell. This Hamiltonian is invariant under the translation of any particle 
coordinate Xk of a length mL, with m integer. This invariance can be expressed as 

77L,ffe(m)]=0 (3) 

where Tk{m) = exp [—j-mL ■ Pk] is the L-translational operator of particle k. Translational invariance implies that 
the eigenfunctions 5* (xi, . . . ,xn) of the cell Hamiltonian eq.© can be written as 

'ffe(TO)*e,_„^ = exp i-iek,m) (4) 

where 9k,m is the eigenvalue associated to the translation Xk Xk — mL. Because of the property of translational 
operators Tk{mi)Tk(m2) = Jfe(mi + 7712), the eigenvalues satisfy 9k,mi + Ok,m2 — 9k.mi+m2- This means that we can 
associate the translational invariance of periodicity L for particle k with an eigenvalue 9k such that 

9k,m = ■m9k (5) 

Let us define an auxiliary wave function as 

$ (xi, . . . , xn) = exp (^^"^-^^k^ * (xi, . . . , xn) (6) 

Using the fact that \l/ is an eigenfunction of the translation operator eq.Q, we get 

^{xi,...,xn) = exp ^-i^{xk - mL)j exp{-im9k)'^ {xi, . . . ,xn) 

= ^{xi,...,Xk-mL,...,XN) (7) 

We have shown that, if 4* is an eigenfunction, then the function $ defined by eq.® is a periodic function of period 
L. The same reasoning can be done for any particle k = 1,...,A^. This means that the eigenfunctions of the 
translationally invariant Hamiltonian eq.(l2]) can be written as 



(xi, . . . jXjv) = exp I i—Y9kXk | $ (xi, . . . , xjv) (8) 

where $ is a periodic function, that is invariant under the translation mL of any particle coordinate. 

Since the Hamiltonian eq.© is invariant under the translation of any particle coordinate separately, one could 
consider in principle a different phase 9k for each particle. However the indistinguishability of nucleons imposes that 
all the phases must be equal, as we now show [4] . Let us consider a phase 9i for the L-translation of particle 1, and a 
phase 92 for the L-translation of particle 2: 

* (.Ti + i,.T2, . . . jXat) = exp{i9i)'^ {xi,X2,...,xn) (9) 
{xi,X2+ L,...,xn) = exp{i92)'^ {xi,X2,.--,xn) (10) 

Applying the permutation symmetry to eq.® gives 



{xi,X2, ■ . ■ , Xn) = - exp {-i9i) * {x2, xi+ L,. . . , xn) 



(11) 
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Applying a translation — L to the second coordinate gives 

(a:i,X2, . . .^xn) = -exp(-z (6*1 - 6*2)) * {x2,xi,. . .,xn) (12) 
Applying the permutation symmetry once again we get 

'^{xi,X2,-.-,Xn) =exp(-j(6'i -e2))^^{Xl,X2,■.■,XN) (13) 

which shows that the two phases must be equal, Oi = 62 = 6l- To be precise we could have 9i — 62 = 2n7r with any 
integer n, but the phase 9l can be taken without any loss of generality in the interval (0, 27r] or (— tt, it], that is in the 
first Brillouin zone. Indeed if we consider a very large number Nr of replicas of the WS cell, 00, the effect of 

the Bloch phase is negligible and we can write global periodic boundary conditions for the wave function 

^ {Xi +NrL, . ..,XN+NrL) = ^ {xi,.. .,xn) (14) 

This implies exp{i0LNr) = 1 or 9l = (2?i — Nr)Tr/Nr with n integer. Let us write n = n'Nr + n" with 1 < rt" < iV^ 
and n' integer, then 

0L = {2n' -1)tt + 6 (15) 

with — TT < 9 < TT. Finally the wave function is 

f 9 ^ \ 

* (a;i, . . . , a; at) = exp I X! * ^^i' ' • ' '^^^ ^^^^ 

This Bloch formulation is a possible way to represent the eigenfunction of the exact N-body Hamiltonian in the cell, 
accounting for the infinite range antisymmetrization correlations with the replicated system. It amounts to introducing 
an extra quantum number 9, 



* (xi, . . . ,a;Ar) = {xi,. . . ,XN\n,9) (17) 

where n denotes the ensemble of the other good quantum numbers. The system observables will thus explicitly 
depend on the phase 9. 

The form of the periodic Hamiltonian eq.(IT]) implies that the wave function of the global system can be expressed 
as an antisymmetrized product of A^-body Bloch wave functions (|16p according to 

00 

*tot (Xl,..., Xoo) - i n *«™«'---'^^) (18) 
m=l 

The Twist Averaged Boundary Condition method consists in assuming that the different angles 9m in the product 
are all different. Then the physical value of any observable O can be computed from eq. ([T8| with the extra restriction 
9m ^ 9n giving 



E^"'^'"!^!".^ (19) 

m— 1 

We can see that within this decoherence hypothesis among the different phases, possible non-diagonal terms disappear 
and the observables per unit cell are obtained as simple averages over phases 

{d)i:^^ r d9{rX9\d\n,9) (20) 
27r 

We will explicitly show in the following that the parts of the Fock space corresponding to different phases are 
indeed disjoint in the case of Slater determinants. For more complex eigenstates. eg . has to be considered as an 
approximation which however appears very well verified in practical applications [3|. 
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B. Application to Slater determinants 

In the following we look for an ansatz for $ (xi, . . . , xat) to be introduced as a variational approximation to the full 
N-body problem. In particular mean field approaches, including Antisymmetrized Molecular Dynamics (AMD) and 
Fermionic Molecular Dynamics (FMD), are based on a Slater approximation. This means that the variational ansatz 
can be written as 

N 

^{xi,...,XN)^A'[[uk{xk) (21) 

fc=i 

where A is the antisymmetrization operator among the N particles and Uk{x + mL) = Uk(x) Vfc. To fulfill this 
periodicity condition we can write 

oo 

Uk{x)^N gk{x~mL) (22) 

m— — oo 

where g is an arbitrary function and A/" is a normalization factor. In the following we will show applications with the 
non-orthogonal single-particle FMD/ AMD basis set[2^[23| given by no n- normalized gaussian wave packets 

gz,{x)=exp(^^^jZk-xf^ (23) 

where Zk, au are complex variational parameters. For this specific choice the translation of the argument is equivalent 
to a translation of the gaussian centroid Zk 



Uk{x) = lim -j= V 9Zk+mL{x) (24) 

^ ^ ■m=-Nr/2 

where the dependence on is implicit and omitted to simplify the notations. It is important to remark that the 
choice of the normalization in eq. ([24l) guarantees that the norm of the wave function is finite, which will allow the 
numerical evaluations below. In the special case where the gaussian width is sufficiently small respect to the cell size, 
such that the overlaps between gaussians can be neglected this norm is readily evaluated as < Uk\uk >= ^y^^\ak\■ 
Let us consider a generic one-body operator O — X^^i "^fc- The matrix element in r-space representation reads 

Nr T 

/OO /"OO f f'^ 

dx I dx'u*{x')o{x' , x)uk{x) = lim / dx dx'u*{x')o{x',x)uk{x) (25) 
-oo J-oo N^^coJ_i^^ 

Because of the periodicity of the Uk this simplifies to 



L/2 
L/2 



Ojk = j^l™ Nr j dx j dx'u*{x')o{x' ,x)uk{x) (26) 



where we have used the fact that the operator is translationally invariant. For a local one-body operator A — X^fcLi '^k 
the matrix element simplifies to 

Gjk = lim Nr / dx u*{x)a{x)uk{x) (27) 

This shows that, because of the periodicity, only integrals over one single box are needed. The expectation value of 
the observable O is given by (for the moment we are ignoring the Bloch phase, and considering expectations only over 
the auxiliary function $, which we indicate <>$): 

N 

< 6 >$=< $|6|$ >= ''kk (28) 
fe=l 

if the Uk constitute an orthonormal basis, and 

j,k=l 
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if not. Here Bjk is the overlap matrix 

Bjk =< Uj\uk > ■ (30) 
Implementing the ansatz ^M^ . the matrix element of a generic local operator is calculated by 

N N^"^oo^WW ^ ^ \ dxgl^^„^^{x)a{x)gz,^^'L{x) (31) 

Because of the finite norm of u/c, the double sum is a convergent quantity with increasing number of replicas 

^^.■/2 K/2 

linijv,,Ar;->oo X! ^ / '^2:g|^+™L(a;)a(a::).gZfc+m'L(2;) = (32) 

Arr/2 w;/2 

=< Ufc|a|Mfc > (34) 
This means that the limits in ea. pil) can be performed separately giving: 

rL/2 
L/2 



limjv.,iv:-^oo — XI XI / '^^3z,+mL(a;)a(a;)gZfc+m'L(a;) (33) 



/L/2 oo 
c^a^ 5z,+mL(2;)a(a;)5z,+m'L(2;) 
-L/2 , 

' m,7n — — oo 



(35) 



In the hypothesis that the spatial extension of the wave function is negligible beyond a linear size L^ax — ML, 
limx^±L^^^{x\gz) = 0, we can write 

^L/2 A/ 

ajk = / dx Y] 9*z,+mLixMx)gz,+,n'L{x) (36) 

' m,m' — — M 

which is a feasible integral. Ea. (|55)) is readily generalized to the case of a non-local operator. 

Because of the translational invariance the wave function has to be multiplied by the Bloch phase factor ea. (IT5|) . 
This means that the expectation values of observables have to be corrected respect to eq. (P5|) . (P^ 

< 6 f: V. - 11^0]^ (37) 

j,k=l 

where the notation Oxj, corresponds to a specific choice for the (arbitrary) Bloch phase — tt < < tt, 

N 

'^0{xi,...,xn) ^AY\_'^ke{xk), (38) 
fc=i 

■>Pke (x) = exp {iOx/ L) Uk (x) (39) 

We can finally write the expression of the different matrix elements for the infinite periodic system, for a given choice 
of the Bloch phase 6. The overlap matrix is not influenced by the Bloch phase, and the same is true for the matrix 
element of any local operator: 



M ^L/2 

"ifce = / dxg*z^^^Lix)a{x)gz^+m'Lix) =ajk (40) 

In particular the local density 

A 

Pe{x) = Yl <x\il^k><'ikj\x> B~J 

l.k=l 
A 

= X < x\uk >< Uj\x > B^^^ = p{x) 

3,k=l 
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is independent of 9. 

The situation is different for non-local operators, because in this case the Bloch phase does not cancel any more 



Ojke ^ X! X! / / <^x'9*z^+7nLix')o{x,x')gz^+„i'L{x)exp{i9{x - x')/L) (41) 



m= — oo m' = — M ' 



Let us take the exemple of the momentum operator 



kiPj{x) ^ '3xp(i6la;/L) f — + — j gz^.+„L(a;) (42) 



The expectation value of the total momentum in the cell K — given by 



AT 



ji=l 

N 9 9 

= J2 < "j l^l"^ > +N-^<k >$ +N- (43) 

We can see that the Bloch phase physically represents a boost to the center of mass motion. It will cancel when 
performing the average between the different phases —it < 9 < tt (TABC), but this will give a finite contribution to 
the kinetic energy. Indeed if we consider the square momentum 

*^ / 2i9 d 9^ \ 

P^,{x)^-expit9x/L) J2 ( 5^ + — (44) 



m=-M 

This gives rise to a total kinetic energy expectation value 



29 ^ ft2 ^2 



h^9 N fh^9^'^ 



This extra kinetic energy term linked to the Bloch phase explicitly enters in the energy variation (or in the equations 
of motion, in the case of dynamical models (2^ [23l . [26l [27j). If we consider the typical size of a Wigner-Seitz cell 
as calculated by Negele and VautherinflJ], in the inner crust {L ^ 15 fm) the phase effect gives rise to an energy 
contribution Ai? « 11 MeV, which is far from being negligible if we compare to the corresponding Fermi energy 
Ep w 24 MeV [iMl. 



C. One-dimensional model cases 



As a first model case, we consider one-dimensional free particles. Even if this system is clearly very far from 
the correlated dishomogeneous solutions relevant for stellar matter, it has the advantage of being exactly solvable. 
Moreover, it is an especially interesting test case for molecular dynamics models [2^ [26l l27j and more general models 
employing localized wave functions[23]: these models are optimized to treat density fluctuations but not necessarily 
adapted to treat the coupling to the continuum which is needed for the extreme neutron-proton ratio which is found 
in the inner crust of neutron stars. 

In one dimension, the energy per particle of a degenerate ideal Fermi gas at density p = N/L = is given by 

epG — ' Qjy^ ■ In the case of a finite system of N particles within a length L with periodic boundary conditions, the 
energy levels are given by 

f2Tm,V , , 
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where is an integer within the interval [—^,^]- Due to this level degeneracy, the total cell momentum (K) is 
equal to for an even number of particles, while it is zero for an odd number. 

Because of the periodicity of the system, it is easier to work in Fourier space. All one-body observables can be 
expressed as a function of the one-body density in momentum space p^(/c,fc') = J2fj=i{^\'4'i){i^jW)B~j^, where the 
single-particle states are given in the momentum representation. These states can be obtained by taking the Fourier 
transform of the wave function. Considering that a periodic function can always be expressed as a Fourier series, 

oo 

u^ix) = c„ze(-2-'-/^) (47) 



1 — — 00 



we get 



^no{k) = V2^ CniS{e/L-k-2nl/L) 

1 — — QO 

= lim V V .f^y dxgz„+mL{x)e^'^'^/'^5{9/L-k~27rl/L) 



oo I WrJ^ 

= lim y \ —- [ ' dxgz ix)e^''''^/^S (0/1 - k - 2Trl/L) 

Nr^OC V Nr L / NrL 

^ ^./^ Y e-'^('^ye''^^-5{e/L-k-2^l/L) (48) 

t — — oo 

At variance with the spatial density, the momentum density thus explicitly depends on the value of the phase 
e = Lko: 

pe{k) = lim Y i?^,py^e-^(^)^'^(^"-^;)<5 {0/L - k - 2^1/ L) (49) 

JVr— s-oo I\j.L ^ — ' 
n,p— 1 

The momentum density is obtained by averaging over the different phase values, giving 

* N 

^*(^) = /5^oo#i E E 54>^e-^(^re^^(-"--;) (50) 

We can see that the value of p in k is solely determined by the wave function corresponding to the phase value 
9 = k — ^ , which corresponds to a unique value in the first Brillouin zone. This means that the number of points 
chosen for the discretization of the integral over the Bloch phase defines the resolution of the momentum density: the 
Bloch phase generates new plane waves which 'fill up' the distribution as needed to obtain a continuous Fermi sea 
starting from the discrete levels of a finite system. 

The total energy and momentum easily follow 

E) = — / dkp^(k)k^, Ik) = / dkpq,(k)k. (51) 
/i- 2m J _^ \ /* J_^ 

The result for the FMD model is shown in FiglT] We have chosen a particle density p = The size L of the 

elementary cell is directly proportional to the number of particles in the cell according to L = The FMD ground 
state is obtained minimizing the total energy with respect to the variational parameters z = {Zn,an,n = 1, . . . , N}. 
This can be obtained using the gradient method, or alternatively by the numerical solution of the FMD equations of 
motion [2^ 
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It is possible to show[26j that if a friction term is added as a multiplicative factor on the r.h.s. of ea. ((52l 

d{H)^ 



dzi 



(54) 



it gives a dynamics of the variational parameters leading to a decrease in time of the total energy. The resulting 
dissipative dynamics can be written as: 



d{H}^ 
dt 



(55) 



and this equation is numerically followed until convergence. 



„^ 
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FIG. 1: FMD calculation of a system of three free particles in a one dimensional box of size L with Bloch boundary conditions. 
Upper part: behavior of the total linear momentum (left) and total energy (right) as a function of the Bloch wave number 
ko — 9/ L Full lines: total expectation values {K)^, (-ff)* (see text); dashed lines: averages over the variational part of the 
wave function {K)is>, {H)<s> without the contribution of the Bloch phase (see text); dash-dotted line: total energy average over 
the Bloch phases. Lower part: spatial (left) and momentum (right) density. 



As expected, Figure [T] shows that the periodicity in fc-space is reproduced by the calculation. This is due to the 
implicit dependence of the expection values (f")*, (H)^ on the choice of the Bloch quantum number fcg = 9/ L obtained 
through the variational procedure. 

The exact results of a free Fermi gas are recovered in the model. This means that the FMD ansatz for the wave 
function appears adapted to describe the free particle system. This was not a-priori obvious because the choice of 
single-particle wave packets represents a restriction of the complete Slater determinant variational space, and the 
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plane wave is obtained only as a mathematical limit a — ^ oo of the wave function, which is not exactly accessible in 
numerical simulations. This implies that the FMD model, which has been introduced essentially to describe cluster 
degrees of freedom, can also address continuum states. 

The other interesting point is that the exact Fermi gas result is obtained with any arbitrary number of particles in 
the simulation. For the 3-particles system shown in Figure [1] the triangles in the lower right part of the figure show 
the momentum values allowed for this finite system. These values are eigenvalues of the periodic part of the wave 
hmction and their spacing decreases with increasing number of particles. As shown by eg. ([50]) . the effect of the 9 
quantum number is to produce all the other missing values for the k quantum number which would be accessible to 
the infinite system, thus reproducing the continuous Fermi distribution. 

In conclusion, the application of TABC to Bloch single particle wave functions appears as a very powerful method to 
recover the correct kinetic energies and wave functions at the thermodynamic limit with variational methods applied 
to small systems even in one dimension. A word of caution is however necessary: to demonstrate the Bloch theorem 
we have explicitly used the fact that the states are eigenvectors of the many-body Hamiltonian, which is in general 
not true when using variational methods. For this reason we will develop an alternative scheme in the next section. 

III. THE REPLICA METHOD 

The twist averaged boundary condition are currently applied in QMC calculations of correlated electron systems 
at the thermodynamic limit. In the specific framework of fermionic molecular dynamics, an alternative method 
to deal with the infinite range antisymmetrisation correlations has been proposed, based on the localized Wannier 
representation ^2^]. Taking advantage of the nested structure of the overlap matrix in a periodic system, an analytical 
method of inversion of this matrix is proposed in ref. [13] . We present here a slightly different derivation of the same 
equations, for an application to any arbitrary single particle basis. 

A. General formulation 

Let us consider that our WS cell contains already many different replicas of the physical system, Ntot = NrN where 
finally we will let — > oo. Similarly, we now interpret the simulation cell Ltot = nL as a huge length, n 3> 1, such 
that working with a finite Lt^t can be considered as equivalent to the thermodynamic limit, meaning that the scaling 
has to be fulfilled for all observables O 

<d>=Nr<6>^ (56) 

with 

-, Ntot 

< O ^ E (57) 

independent of the periodicity properties of the Hamiltonian. If we consider the Afot-body wave function, the invari- 
ance property of the system under the translation mLtot of any particle coordinate is still verified. The difference 
respect to the previous section is that the system composed of Ntot particles has an additional symmetry which is 
not directly linked with the use of a finite box L but is rather due to the physical periodicity of the pasta structures 
over a lenght scale L — Ltot/n. We have already observed that, because of the absence of an external field, the cell 
Hamiltonian eq.([2|), which now represents the full Hamiltonian eq.([T|), 

Ntot oo Ntot 

H = '^ii+ ^ '^v{xi- Xj - mLtot) (58) 

i—l 7n— — Qo iyj 

is invariant under the simultaneous translation of every particle coordinate of any arbitrary length. The translational 
invariance which will be relevant to us is the invariance under the simultaneous translation of every particle coordinate 
of the physical length L, defined as the periodic length of the one body density 

p{x + mL) = p[x) (59) 

If the Hartree-Fock field as an external field, this translational invariance would be the only one verified by the mean- 
field hamiltonian, and it would be verified for each particle coordinate separately. The situation is different in the 
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case of an interacting system as described by the self-consistent Hartree-Fock approach. In this case, translational 
invariance with respect to any arbirary length is respected because of the self-consistency of the mean-field approach, 
Uij = SEhf I^Pji where Ehf is the interaction part of the Hartree-Fock energy. However, this translational invariance 
applies only to the simultaneous translation of all particle coordinates. The only special feature of the length L is 
that eq. (|59p is fulfilled, which will allow us to impose a simplified expression for the one-body wave functions, as we 
will see in a moment. 

Following the derivation of the previous section, these two invariance properties lead to the definition of a Bloch 
phase for the TV^ot-body wave function according to 

* (xi, . . -.XNt^t) = exp X! ■ • ■ '^^tot) (60) 

\ *°* fc=i / 

where $ is a periodic function, that is invariant both under the translation mL of all Nt^t particle coordinates 
simultaneously, and under the translation mLtot of any particle coordinate. 

The difference respect to the derivation of the previous section is that now we have a factor 1 / Ntot multiplying the 
Bloch phase. Since we are at the thermodynamic limit Ntot — >■ oo, the correction on kinetic observables induced by 
this Bloch phase can be neglected and we can consider that the total A^tot-body wave function shares the periodicity 
properties of $. 

Let us now look for an ansatz for this L— and Ljot— periodic \l/ (xi, . . . ,XNtot) to be introduced as a variational 
approximation to the full A^'tot-body problem. We take the same Slater ansatz used in the previous section: 

Ntot 

<i/{xi,...,XN,„t) = AY[uk{x) (61) 

fe=i 

with Uk{x + mLtot) — Ukix). Because of the periodicity of the one body density, each sub-cell of linear dimension L 
can be exactly associated to a finite number N — Ntot/Nr particles, and the wave function can be written introducing 
a sub-cell index m as: 

N 

*(xi,...,xjv.„j =in n ) (62) 

k—l m — 1 

Moreover, the periodicity of p can be written as 

Nr N 

p{x + L) = ul„^{x + L)uj^„,.{x + L)BJ^ 

m,m' — 1 k.j — 1 
Nr N 

= Yl ^ Uk,m(.x)Uj,m'{x)B~,]- = p{x) (63) 
m,m' — l k.j — 1 

This condition can be satisfied introducing only N independent one body wave functions Uk^m{x) = U}.{x — mL) where 
in principle Uk are arbitrary (non periodic) functions. If these functions would form an orthonormal basis 

< j, m|fc, m' >^ SjkSmm' (64) 



the problem would be reduced to a body problem as with the ansatz eg . (f2T|) . Indeed in the calculation ([571) of an 
arbitrary one-body observable only the Uk, k = 1, . . . , N functions would be needed 

^ Ntot N 

<d>^= — Yokk = Y^kk. (65) 

k=i fc=i 

B. Application to non-orthogonal single-particle basis 

For practical applications we will not work with an orthonormal basis, but we want to take the FMD/AMD gaussian 
wave packet variational ansatz eg. (|23|) . namely 



Uk{x - ml) = gZk+ral{x) 



(66) 
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Since the gz^+mi do not form an orthonormal basis, we are left with a problem of computational size N^^^, that is 
impossible to solve (recall Ntot — >■ oo). 

A possibility would be again to impose the simplified expression: 



^{xi,...,XN,J = n ■^Haz.+miix) = n '^{xi,...,Xn) (67) 

m—1 k—1 m—1 

which allows to work with only N functions. 

The conceptual problem with this simple ansatz eg . ([ST]) is that it neglects the antisymmetrization correlations 
among different cells. We know that such correlations are important as they are at the origin of the band structure. 
For this reason we need to keep the antisymmetrization correlation among all the Ntot = Nr ■ N particles. 

By introducing a linear transformation 

uk = Ugl (68) 

defined by 

Uk,n{x) = ^= expUm'^^ gz^+mL{x), (69) 



we can reduce the computation size to N^N'^ keeping the antisymmetrization correlations among all the Ntot — NNr 
particles. It is interesting to remark that, in the limit of iV^ — > oo, this transformation coincides with the definition 
of Wannier functions [l']. This transformation preserves the determinant 

^ n 9Z,+mLix) = i n (70) 
m—1 n—1 

so that we can write 

Nr- N 

'ifixi,...,XNtJ^AYlYluk,n{x) (71) 



n=l k=l 



The advantage of the transformation l|69|) is that matrices are block-diagonal, that is functions corresponding to 
different phases 7 = (2n — Nr)7r/Nr are orthogonal 

< Uj^n\uk,n' >= 5n.n'Bjk{n)- (72) 



We have thus shown that the decoherence hypothesis eq. (|T9)) of the TABC method is exactly verified in the case of 
variational approaches based on Slater determinants. 

Having recovered the orthogonality on the level of the replica quantum number we can write 

^ Nr N 

<0>%= ^^Y.Y. < ^i.n\o\uk,n > Bfc/(n) (73) 

n=l j,k=l 

The matrix element of any operator is given by 

/L/2 oo / 2nn\ 

dx 9*z,+mLix')o{x,x')gz^+m'Lix)exp(-i{m-m')j-^\ (74) 

and explicitly depends on the n quantum number associated to the antisymmetrization among the different subcells. 
In particular the overlap matrix reads 

dx Y 9*z,+mLix)9Zk+m'Lix)expi-i{m~m')j-^] (75) 

Using the same argument that we had for the expectation (|36p . the series can be reduced to a finite sum 

, 27m 



M .L/2 / 

Bjk{n)= Y / dxg*z^+„-,L{x)gz^+m'L{.x)exp{-i{m-m') 

m..rn.' = -M '' '^/"^ ^ 



m,m' — — M ^ 



LNr 



(76) 
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In practical application, a large value of M « 10 is needed to describe non-interacting particles, but M = 1 turns 
out to be sufficient in the presence of the nuclear and Coulomb interaction for the particle densities relevant for the 
neutron star crust. From the computational viewpoint, M = 1 is equivalent to calculate a single gaussian overlap 
between each pair if periodic boundary conditions are applied to the gaussian, meaning that the simulation of the 
infinite system is not more expensive numerically than the one for the finite system. The analogous expression for 
the matrix elements reads 



< Uj,n\a\Uk^n > = 



M 

E 



m,m' — — Af 

Going to the Nr oo limit we can write 



with — TT < 7 < TT. The observables read 



LNr 



L/2 / 2j^ji 

dxg*z^+mLi^)a-{x)gZk+,n'L{x) exp ( -i{m - m') 

-L/2 \ 



27rn 
hm 

AT^^OO LNr 



< o >i= 



' J — TT -I 1 



(77) 



(78) 



(79) 



which is expected to give equivalent results to the implementation where TABC are applied. 




FIG. 2: Energy associated to a system of two free FMD particles in one dimension at a density p = as a function of the 
number of replicas. Dash-dotted line: free Fermi gas in one dimension. 



The FMD ground state energy for a system of two free particles in one dimension is shown in Figure [2] as a function 
of the number of replicas, with the Wannier choice for the wave packets. For each calculation, the number M was 
chosen high enough to have convergent results. We have already observed that the FMD variational space is sufficiently 
flexible to describe a free-particles system. The FMD solution for each value of Nr thus represents the exact energy 
of a finite system with periodicity L — 2/'p and size Ltot — L ■ Nr- The thermodynamic limit, corresponding to the 
anlytical Fermi gas result, is obtained for about ten replicas, corresponding to ten values for the phase 7 of eg. (fTSt . 



IV. COMPARISON OF THE TWO METHODS 



In the previous sections, we have introduced the Bloch and the Wannier representation as two alternative ways to 
address the thermodynamical limit in variational theories. From a principle point of view the two approaches are very 
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similar, however they lead to very different expressions in the calculation of observables. In particular, the fact of 
considering only diagonal terms in the phase quantum number naturally emerges as a consequence of the orthogonality 
of the Wannier basis in the replica method, while it appears as a simple working hypothesis in TABC. Moreover, 
the equivalence is expected in the ideal case of an infinite number of twists and under the condition that the wave 
function is an exact eigenvalue of the Hamiltonian. This first condition is numerically expensive, while the second is 
not assured in variational methods. For all these reasons, the equivalence of the two formalisms is not a-priori clear 
and will be discussed in the present section. 

A. Finite number of replicas and finite number of twists 

We now show that, in the hypothesis that the variational theory produces the exact eigenvectors of the many-body 
Hamiltonian, the two methods to treat the antisymmetrisation correlations can be mapped on each other for any 
system size. 

In the Bloch formalism, if we impose that the condition of periodic boundaries ea. (jl4p is verified for a finite 
translation of size N^L, the original problem of an infinite system with periodicity L is transformed into the problem 
of a finite system of size N^L, with periodic boundary conditions. 

The allowed values for 6 are then discrete 9 = ^""^'' tt, where n is an integer. This amounts to discretize the 

integral (PO]) with /S.9 = j^. We expect then the two methods to give exactly the same results. This expectation is 
confirmed by Figure[3J which compares for the case of the non-interacting one-dimensional system the energy obtained 
considering a periodic system of Ntot — 'pNrL particles, with the calculation of a reduced N = system averaged 
over Nr phases with the two (Bloch and replica) proposed methods. We can see that the agreement is perfect for 
all sizes, and all calculations converge to the free Fermi gas result with increasing number of replicas (respectively: 
twists). 

Ltot[fm] 

2.3 6.9 11.5 16.1 20.7 25.3 29.9 34.5 39.1 43.7 48.3 52.9 57.5 

15 
14.5 

14 
13.5 

CD 

Ella 

12.5 

12 
11.5 

11 

1 3 5 7 9 11 13 15 17 19 21 23 25 

Nr 

FIG. 3; Energy of a free periodic one-dimensional system at a density p = as a function of the number of replicas. The 
upper abscissa represents the total linear size in the case of simple periodic boundary conditions (full line), the lower one is a 
number of replicas for the replica method (crosses), and a number of twists for the Bloch method (circles). Dash-dotted line: 
free Fermi gas in one dimension. 

This shows the equivalence of the different methods when the wave vector is an eigenvalue of the Hamiltonian. 
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B. Limitations of the Bloch method in variational applications 



In realistic applications of interacting fermion system, the variational solution is always an approximation of the 
exact eigenvalue. This may create a bias in the different method that we turn to discuss. For eq.Q to correctly 
account for the correlations due to the particles outside the cell L, it is necessary that the wave function Vl/g is an 
eigenfunction of the hamiltonian. Conversely the Wannier representation can be obtained from the solution of the 
infinite periodic system via a simple linear transformation eq. (|68p . Of course if the variational ansatz is not adequate 
to describe the exact state both representations will be false. However the error due to the inadequacy of the solution 
in the cell L will be propagated in an uncontrolled way in the Bloch method for the description of the global system, 
and this will not be the case for the replica method that directly addresses the replicated system of size TV^L. 

To illustrate this difference, we can once again take the simple example of the non-interacting Fermi gas. There is 
a popular microscopic theory in nuclear physics, the so-called AMD modeljl^l, where the single particle wave packets 
are gaussians of fixed width. With this ansatz the variational space does not contain the exact solution of this problem 
at the bulk limit, that is an antisymmetrized combination of plane waves. In particular, if we choose a small value for 
the width, as it is variationally obtained to describe finite nucleifl^l, it is clear that the states will be very far from 
the exact eigenvectors of the kinetic hamiltonian. The comparison of the two methods in this model case is shown in 
Figure IH 



20 



15 



10 




FIG. 4: Comparison of the different boundary conditions for the AMD model in the case of a non interacting system as a 
function of the system size. Dash-dotted line: free Fermi gas. Dashed line: AMD ground state energy for a system size large 
enough for finite size effects be negligible (L = 150). Circles: finite system with Bloch boundary conditions. Cross: replica 
method. Stars: periodic boundary conditions. 

We can see in Figur^ that the AMD ansatz is not adapted for this problem, from the fact that asymptotically 
the energy of the free Fermi gas (full line) is not recovered by the AMD calculation (dashed line). This asymptotic 
value is perfectly reproduced by the replica method (crosses) for any system size, while the Bloch method (open 
circles) produces artificial fiuctuations and attains convergence only when the Bloch phase is negligibly small, and the 
calculation is perfectly equivalent to simple periodic boundary conditions (stars). 

The non-interacting system is certainly not the ideal application ground of molecular dynamics model. As a second 
model case which is closer to the physical case of the neutron star crust we take a one-body periodic potential, which 
is obtained as a periodic self-consistent mean field in realistic applications of mean-field variational models [THlOl [l3 | if 
the two-body nuclear and Coulomb interaction is added, including the interaction with a uniform electron background. 
For this model one-dimensional application we only wish to discuss the effect of the boundary conditions, therefore 
we restrict to the simpler case of an external periodic potential. In particular, the choice of a gaussian potential has 
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the advantage of making ah calculations analytical. Let us take a potential form: 



Vgaus [x) = VJ) ^ ^ exp 

^■. — _N_ m — — OO 



Za 



(80) 



where Vb, a, &„ are parameters. In particular, the choice a = I, 6„ = ^ for a system with average density p, corresponds 
to having a single particle per potential well. The average energy is readily calculated as 

N 

. \ I 

^gaus 



1 ^ / 

Vo lim — y Br\ bn 



hJ=l 



a*a + aiU + a^a* 



E 



exp 



mi ,m2 ,n— — oo 



a*{Zi - bn + niiLy + ai{Z* - bn + ^2^)^ + aiZt - Z* - {mi - m2)L) 

a*a + a^a + 0^0* 



Results from the Bloch method for a case of 3 particles with an average density p = are presented in Figure [S] 
At variance with the free particle case, the periodic part of the wave function $ now explicitly depend on the 9 
quantum number, leading to a dependence of the one body density. 

The comparison between the different methods of treating the boundary conditions is shown in figure [51 
This model is less trivial than the free Fermi gas, although still very schematic. The FMD model is expected to 
be a good approximation of this system, but there is no guarantee the FMD solution should be exact. Similar to the 
previous application, the replica method gives by construction the asymptotic result for any finite size, provided a 
sufficient number of phases is considered. Conversely for the Bloch method, which is calculated with the same number 
of phases as the replica method for the application of FiglHl a convergence is reached only when the effect of the phase 
can be neglected. 

This model example shows that in realistic applications, where any variational ansatz might be simply an approxi- 
mation of the exact energetics of the system, the different boundary conditions are not equivalent, and the employ of 
Bloch corrections to the single particle basis may increase the deviation respect to the exact solution. 



V. CONCLUSIONS 



In this paper we have introduced and critically discussed two different ways of addressing the boundary conditions 
in finite fermionic system in order to account for the infinite range antisymmetrisation correlations present at the 
thermodynamical limit. In the case of state vectors which are eigenstates of the many-body Hamiltonian, we have 
shown that the use of Bloch wave functions with twist averaged boundary conditions physically corresponds exactly to 
a bigger system constituted of a number of replicas of the system under study equal to the number of phases. In turn, 
this last formalism can be interpreted as the use of localized Wannier states. While these different representations are 
equivalent when dealing with eigenvectors, the same is not true in the case of approximate solutions of the many-body 
problem, as it is the case in the variational approach. In this case, the approach of the thermodynamic limit is not 
properly treated by the Bloch method, while the replica method provides a very accurate evaluation of the Fock 
energy. Even in the case of a non-orthogonal single particle basis this method can be numerically implemented with 
a moderate computational cost, opening the possibility of a completely quantum-mechanical treatment of nuclear 
matter with molecular dynamics approaches. 
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